# This code replicates all analysis of survey data including Tables 1, 3, A2, A6-A9 and Figures 3, 5, A4, A6-A11, and A15-A17

library(estimatr)
library(dplyr)
library(estimatr)
library(texreg)
library(glmnet)
library(ggplot2)
library(patchwork)
library(upstartr)
library(RColorBrewer)
library(magrittr)

set.seed(416)

load("data/survey_data.Rdata")

# Table 1 -----------------------------------------------------------------

get_dif <- function(model){
  e <- coef(model)["Z_common"] - coef(model)["Z_alt"]
  se <- sqrt(vcov(model)["Z_common", "Z_common"] + vcov(model)["Z_alt", "Z_alt"] -
               2 * vcov(model)["Z_common", "Z_alt"])
  return(2 * pt(-abs(e/se), df = model$df["Z_common"]))
}


m1 <- lm_robust(past_meeting ~ Z_common + Z_alt + as.factor(block_ID), 
                data = survey, clusters = nro_cuadra)
m2 <- lm_robust(past_meeting ~ Z_common + Z_alt + as.factor(block_ID), 
                data = survey, clusters = nro_cuadra, subset = in_endline == 1)
m3 <- lm_robust(attended_past ~ Z_common + Z_alt + as.factor(block_ID), 
                data = survey, clusters = nro_cuadra)
m4 <- lm_robust(attended_past ~ Z_common + Z_alt + as.factor(block_ID), 
                data = survey, clusters = nro_cuadra, subset = in_endline == 1)
m5 <- lm_robust(any_exposure ~ Z_common + Z_alt + as.factor(block_ID), 
                data = outcomes, clusters = CUADRANTE)
m6 <- lm_robust(any_exposure ~ Z_common + Z_alt + as.factor(block_ID), 
                data = outcomes, clusters = CUADRANTE, subset = panel == 1)
m7 <- lm_robust(attended ~ Z_common + Z_alt + as.factor(block_ID), 
                data = outcomes, clusters = CUADRANTE)
m8 <- lm_robust(attended ~ Z_common + Z_alt + as.factor(block_ID), 
                data = outcomes, clusters = CUADRANTE, subset = panel == 1)

ctrl_means <- c(summary(survey$past_meeting[survey$Z_alt == 0 & survey$Z_alt == 0])[4],
                summary(survey$past_meeting[survey$Z_alt == 0 & survey$Z_alt == 0 & survey$in_endline == 1])[4],
                summary(survey$attended_past[survey$Z_alt == 0 & survey$Z_alt == 0])[4],
                summary(survey$attended_past[survey$Z_alt == 0 & survey$Z_alt == 0 & survey$in_endline == 1])[4],
                summary(outcomes$any_exposure[outcomes$Z_alt == 0 & outcomes$Z_alt == 0])[4],
                summary(outcomes$any_exposure[outcomes$Z_alt == 0 & outcomes$Z_alt == 0 & outcomes$panel == 1])[4],
                summary(outcomes$attended[outcomes$Z_alt == 0 & outcomes$Z_alt == 0])[4],
                summary(outcomes$attended[outcomes$Z_alt == 0 & outcomes$Z_alt == 0& outcomes$panel == 1])[4])

ctrl_means <- as.character(round(ctrl_means, 3))

write(texreg(list(m1, m2, m3, m4, m5, m6, m7, m8), omit.coef = c("block|Intercept|Z_alt"), include.ci = F, digits = 3,
             custom.coef.names = c("Meetings"),
             custom.model.names = paste("(",as.character(c(1:8)), ")", sep = ""), 
             custom.header = list("Baseline exposure" = 1:2, "Baseline attendance" = 3:4, "Endline exposure" = 5:6, "Endline Attendance" = 7:8),
             custom.gof.rows = list("Block FE" = c(rep("$\\checkmark$", 8)),
                                    "Observations" = lapply(list(m1, m2, m3, m4, m5, m6, m7, m8), function(x){summary(x)$nobs}),
                                    "Clusters" = lapply(list(m1, m2, m3, m4, m5, m6, m7, m8), function(x){x$nclusters}),
                                    "Control mean"= ctrl_means,
                                    "Sample"= list("Baseline", "Panel", "Baseline", "Panel",
                                                   "Endline", "Panel", "Endine", "Panel"))),
      file = "results/Table_1.tex")




# Table 3 -----------------------------------------------------------------


panel <- filter(outcomes, REGISTRO_BD %in% survey$REGISTRO_BD)
m1a <- lm_robust(I(trust_police - trust_po_bl) ~ any_exposure, 
                 data = panel, subset = Z_common == 1, cluster = CUADRANTE)
m2a <- lm_robust(I(trust_police - trust_po_bl) ~ attended, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m3a <- lm_robust(I(trust_police - trust_po_bl) ~ any_exposure + attended, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m4a <- lm_robust(I(po_index - po_index_bl) ~ any_exposure, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m5a <- lm_robust(I(po_index - po_index_bl) ~ attended, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m6a <- lm_robust(I(po_index - po_index_bl) ~ any_exposure + attended, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)

m1b <- lm_robust(I(trust_police - trust_po_bl) ~ any_exposure + trust_po_bl, 
                 data = panel, subset = Z_common == 1, cluster = CUADRANTE)
m2b <- lm_robust(I(trust_police - trust_po_bl) ~ attended + trust_po_bl, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m3b <- lm_robust(I(trust_police - trust_po_bl) ~ any_exposure + attended + trust_po_bl, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m4b <- lm_robust(I(po_index - po_index_bl) ~ any_exposure + po_index_bl, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m5b <- lm_robust(I(po_index - po_index_bl) ~ attended + po_index_bl, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)
m6b <- lm_robust(I(po_index - po_index_bl) ~ any_exposure + attended + po_index_bl, 
                 data = panel, subset = Z_common == 1,  cluster = CUADRANTE)

# calculate means and controls for unexposed respondents (relative to each column) for footer
means <- with(panel, list(round(mean(trust_police[any_exposure == 0] - trust_po_bl[any_exposure == 0], na.rm = T), 3),
                       round(mean(trust_police[attended == 0] - trust_po_bl[attended == 0], na.rm = T),3),
                       round(mean(trust_police[any_exposure == 0] - trust_po_bl[any_exposure == 0], na.rm = T),3),
                       round(mean(po_index[any_exposure == 0] - po_index_bl[any_exposure == 0], na.rm = T),3),
                       round(mean(po_index[attended == 0] - po_index_bl[attended == 0], na.rm = T),3),
                       round(mean(po_index[any_exposure == 0] - po_index_bl[any_exposure == 0], na.rm = T), 3)))
              
sds <- with(panel, list(round(sd(trust_police[any_exposure == 0] - trust_po_bl[any_exposure == 0], na.rm = T), 3),
                       round(sd(trust_police[attended == 0] - trust_po_bl[attended == 0], na.rm = T), 3),
                       round(sd(trust_police[any_exposure == 0] - trust_po_bl[any_exposure == 0], na.rm = T), 3),
                       round(sd(po_index[any_exposure == 0] - po_index_bl[any_exposure == 0], na.rm = T), 3),
                       round(sd(po_index[attended == 0] - po_index_bl[attended == 0], na.rm = T), 3), 
                       round(sd(po_index[any_exposure == 0] - po_index_bl[any_exposure == 0], na.rm = T), 3)))

# calculate p-values on differences in coefficients across the panels
compare_ests <- function(m1, m2){
  ncoef = length(coef(m1))
  dif = coef(m1)[2:ncoef]-coef(m2)[2:ncoef]
  se = sapply(2:ncoef, function(x){sqrt(vcov(m1)[x,x] + vcov(m2)[x,x])})
  return(round(2 * pt(dif/se, summary(m2)$df), 3))
}


p_exposure <- c(compare_ests(m1a, m1b), "", compare_ests(m3a, m3b)[1], 
                   compare_ests(m4a, m4b), "", compare_ests(m6a, m6b)[1])
p_attend <- c("", compare_ests(m2a, m2b), compare_ests(m3a, m3b)[2], 
                  "", compare_ests(m4a, m4b), compare_ests(m6a, m6b)[2])

# Note that the two panels are compiled by hand in the manuscript. all estimates appear in
# the texreg output.
texreg(list(m1a, m2a, m3a, m4a, m5a, m6a), omit.coef = c("po_index_bl|trust_po_bl|(Intercept)"), include.ci = F, digits = 3,
       custom.coef.names = c("Any exposure", "Attended meeting"),
       custom.model.names = paste("(",as.character(c(1:6)), ")", sep = ""), 
       custom.header = list("Change in trust in police" = 1:3, "Change in police quality index" = 4:6),
       custom.gof.rows = list("Block FE" = c(rep("$\\checkmark$", 6)),
                              "Observations" = lapply(list(m1a, m2a, m3a, m4a, m5a, m6a), function(x){summary(x)$nobs}),
                              "Clusters" = lapply(list(m1a, m2a, m3a, m4a, m5a, m6a), function(x){x$nclusters}),
                              "Untreated mean"= means,
                              "Untreated standard deviation"= sds),
       file = "results/Table_3a.tex")

texreg(list(m1b, m2b, m3b, m4b, m5b, m6b), omit.coef = c("po_index_bl|trust_po_bl|(Intercept)"), include.ci = F, digits = 3,
       custom.coef.names = c("Any exposure", "Attended meeting"),
       custom.model.names = paste("(",as.character(c(1:6)), ")", sep = ""),
       custom.header = list("Change in trust in police" = 1:3, "Change in police quality index" = 4:6),
       custom.gof.rows = list("Observations" = lapply(list(m1a, m2a, m3a, m4a, m5a, m6a), function(x){summary(x)$nobs}),
                              "Clusters" = lapply(list(m1a, m2a, m3a, m4a, m5a, m6a), function(x){x$nclusters}),
                              "Untreated mean"= means,
                              "Untreated standard deviation"= sds,
                              "Any exposure: Panel A = Panel B, p-value"= p_exposure,
                              "Attended mtg: Panel A = Panel B, p-value"= p_attend),
       file = "results/Table_3b.tex")



# Figure 3 (and Tables A6-A8) ---------------------------------------------
## note that Figure 3 and Tables A6-A8 report the same ITT estimates. this code
## generates both Figure 3 and the tables.

## function generates three regression specifications for each outcome
run_models <- function(dv, lagged_dv_string){
  # no covariate adjustment
  formula1 <- paste0(dv, "~ Z_common + Z_alt")
  # includes block FE
  formula2 <- paste0(formula1, " + as.factor(block_ID)")
  # includes block FE and baseline level
  formula3 <- paste0(formula2, " + ", lagged_dv_string)
  
  mod1 <- lm_robust(as.formula(formula1), data = outcomes, clusters = CUADRANTE)
  mod2 <- lm_robust(as.formula(formula2), data = outcomes, clusters = CUADRANTE)
  mod3 <- lm_robust(as.formula(formula3), data = outcomes, clusters = CUADRANTE)
  l <- list(mod1, mod2, mod3)
  return(l)
}

## function calculates mean and standard deviation of an outcome for the control group
get_dv_summary <- function(dv){
  formula <- paste0(dv, " ~ 1")
  m <- lm(as.formula(formula), data = outcomes, subset = Z_common == 0 & Z_alt == 0, na.action = "na.omit")
  mean <- rep(round(coef(m), 3), 3)
  sd <- rep(round(sd(unlist(model.frame.default(m))), 3), 3)
  return(data.frame(mu = mean, sd = sd))
}

## function prepares summaries of dependent variable for table footer
summarize_control_dv <- function(l){
  return(
    do.call("rbind", 
            lapply(l, get_dv_summary)))
}

## function extracts ITT estimates from each model
get_est_df <- function(mod_list, category, outcome){
  df <- do.call("rbind", lapply(mod_list, FUN = function(x) {data.frame(est = x$coefficients[2:3],
                                                                        se = x$std.error[2:3])}))
  return(df %>% mutate(
    treatment = rep(c("Meetings", "Flyers"), 3),
    outcome = outcome,
    category = category,
    model = rep(c("None", "Block FE", "Block FE, lagged DV"), each = 2)
  ))
  
}

## function formats tables
table_maker <- function(lms, ctrl, dv_labs){
  hdr = lapply(1:(length(lms)/3), function(x){(x-1)*3 + 1:3})
  names(hdr) <- c(dv_labs)
  table <- texreg(lms, omit.coef = c("block|_bl|Intercept"), include.ci = F, digits = 3,
                  custom.gof.rows = list("Block FE" = c(rep(c("", "$\\checkmark$", "$\\checkmark$"), length(lms)/3)),
                                         "Baseline outcome" = c(rep(c("", "", "$\\checkmark$"), length(lms)/3)),
                                         "Control mean" = c(ctrl[,1]),
                                         "Control std. dev." = c(ctrl[,2])), 
                  custom.coef.names = c("Meetings", "Flyers"),
                  custom.model.names = paste("(",as.character(c(1:length(lms))), ")", sep = ""),
                  custom.header = hdr, table = F)
  return(table)
}


## Table A6
### Trust outcomes
t1 <- run_models(dv = "trust_police",
                 lagged_dv_string = "trust_bl_imp + trust_bl_miss")
t2 <- run_models(dv = "trust_beat_po",
                 lagged_dv_string = "trust_bl_imp + trust_bl_miss")
t3 <- run_models(dv = "I(trust_beat_po - trust_police)",
                 lagged_dv_string = "trust_bl_imp + trust_bl_miss")
t4 <- run_models(dv = "I(trust_police - trust_po_bl)",
                 lagged_dv_string = "trust_bl_imp + trust_bl_miss")

ctrl_trust <- summarize_control_dv(list("trust_police", 
                                        "trust_beat_po",
                                        "I(trust_beat_po - trust_police)",
                                        "I(trust_police - trust_po_bl)"))
write(table_maker(lms = c(t1, t2, t3, t4), ctrl = ctrl_trust, dv_labs = c("Police", "Beat Police", "Beat Police - Police", "Police - Baseline")),
      file = "results/Table_A6.tex")

### for left panel of Figure 3
trust_df <- rbind(
  get_est_df(t1, category = "Trust",outcome = "Trust in police\n(general)"),
  get_est_df(t2, category = "Trust",outcome = "Trust in beat police"),
  get_est_df(t3, category = "Trust",outcome = "Trust in beat police -\nTrust in police (general)"),
  get_est_df(t4, category = "Trust", outcome = "Change in trust in police\n(endline - baseline)"))

## Table A7
### Police quality outcomes

u1 <- run_models(dv = "po_index",
                 lagged_dv_string = "po_index_bl_imp + po_index_bl_miss")
u2 <- run_models(dv = "po_beat_index",
                 lagged_dv_string = "po_index_bl_imp+ po_index_bl_miss")
u3 <- run_models(dv = "I(po_beat_index - po_index)",
                 lagged_dv_string = "po_index_bl_imp+ po_index_bl_miss")
u4 <- run_models(dv = "fisc_index",
                 lagged_dv_string = "fisc_index_bl_imp+ fisc_index_bl_miss")

u5 <- run_models(dv = "I(po_index_red - fisc_index)",
                 lagged_dv_string = "po_index_bl_imp+ po_index_bl_miss + fisc_index_bl_imp+ fisc_index_bl_miss")

ctrl_qual <- summarize_control_dv(list("po_index", 
                                       "po_beat_index",
                                       "I(po_beat_index - po_index)",
                                       "fisc_index",
                                       "I(po_index_red - fisc_index)"))

write(table_maker(lms = c(u1, u2, u3, u4, u5), ctrl = ctrl_qual, dv_labs = c("Police", "Beat Police", "Beat Police - Police", "Prosecutors", "Police - Prosecutors")),
      file = "results/Table_A7.tex")

### for center panel of Figure 3
quality_df <- rbind(
  get_est_df(u1, category = "Quality",outcome = "(General) police quality"),
  get_est_df(u2, category = "Quality",outcome = "Beat police quality"),
  get_est_df(u3, category = "Quality",outcome = "Beat police quality-\n(general) police quality"),
  get_est_df(u4, category = "Quality", outcome = "Prosecutor quality"),
  get_est_df(u5, category = "Quality", outcome = "Prosecutor quality -\n(general) police quality"))

## Table A8

v1 <- run_models(dv = "conviv_index",
                 lagged_dv_string = "conviv_index_bl_imp + conviv_index_bl_miss")
v2 <- run_models(dv = "I(conviv_index - conviv_index_bl)",
                 lagged_dv_string = "conviv_index_bl_imp + conviv_index_bl_miss")
v3 <- run_models(dv = "sec_index", 
                 lagged_dv_string = "sec_index_bl_imp + sec_index_bl_miss")
v4 <- run_models(dv = "sec_index - sec_index_bl", 
                 lagged_dv_string = "sec_index_bl_imp + sec_index_bl_miss")

ctrl_outputs <- summarize_control_dv(list("conviv_index",
                                          "I(conviv_index - conviv_index_bl)",
                                          "sec_index",
                                          "I(sec_index - sec_index_bl)"))

write(table_maker(lms = c(v1, v2, v3, v4), ctrl = ctrl_outputs, dv_labs = c("Convivencia", "Convivencia - Baseline Convivencia", "Security", "Security - Baseline Security")),
      file = "results/Table_A8.tex")

### for right panel of Figure 3
output_df <- rbind(
  get_est_df(t1, category = "Perceived outcomes",outcome = "Convivencia"),
  get_est_df(t2, category = "Perceived outcomes",outcome = "Change in convivencia\n(endline - baseline)"),
  get_est_df(t3, category = "Perceived outcomes",outcome = "Security"),
  get_est_df(t4, category = "Perceived outcomes", outcome = "Change in security\n(endline - baseline)"))


## generate Figure 3

df <- rbind(trust_df, quality_df, output_df) %>%
  as.data.frame() %>%
  mutate(category = factor(category, levels = unique(category)),
         outcome = factor(outcome, levels = unique(outcome)),
         treatment = factor(treatment, levels = unique(treatment)),
         model = factor(model, levels = unique(model)),
         treatment2 = ifelse(treatment == "Meetings", "ITT of meetings on...",
                             "ITT of flyers on..."),
         treatment2 =factor(treatment2, levels = unique(treatment2)))  %>%
  filter(treatment == "Meetings")


ggsave(ggplot(df, aes(x = outcome, 
                      y =est, col = model, pch = model)) + 
         scale_x_reordered("") + 
         scale_y_continuous("ITT Estimate") + 
         facet_grid(.~category, scales = "free_x") + 
         geom_hline(lty = 3, yintercept = 0) + 
         geom_point(size = 2, position = position_dodge(width=0.5)) + 
         geom_errorbar(aes(ymin = est - 1.96 * se,
                           ymax = est + 1.96 * se), 
                       width = .4, position = position_dodge(width=0.5)) + 
         theme_minimal() + 
         scale_color_manual("Covariates", values = c("#41b6c4", "#1d91c0", "#225ea8")) + 
         scale_shape_manual("Covariates", values = c(15:17)) + 
         theme(legend.position = "bottom",
               strip.text.x = element_text(face = "bold"),
               panel.border = element_rect(color = "black", fill = NA, size = .5),
               axis.text.x = element_text(angle =45, vjust = 1, hjust=1)), 
       file = "results/Figure_3.pdf",
       width = 9, height = 5.5)

# Figure 5 ----------------------------------------------------------------

m1 <- lm_robust(any_exposure ~ -1 + as.factor(trust_po_bl), data = outcomes, cluster = CUADRANTE)
m2 <- lm_robust(attended ~ -1 + as.factor(trust_po_bl), data = outcomes, cluster = CUADRANTE)

fig_5 <- rbind(summary(m1)$coef[,1:2],
               summary(m2)$coef[,1:2]) %>%
  as.data.frame() %>%
  mutate(outcome = rep(c("Any Exposure to Intervention", "Attended Meeting"), each = 4),
         bl_trust = rep(1:4, 2)) %>%
  ggplot(aes(x = bl_trust, y = Estimate)) + 
  geom_point(size = 3) + 
  facet_grid(~outcome) + 
  geom_errorbar(aes(ymin = Estimate - 1.96 * `Std. Error`,
                    ymax = Estimate + 1.96 * `Std. Error`), width = .1) + 
  scale_y_continuous("Proportion of respondents in panel") + 
  xlab("Baseline Trust") + 
  ggtitle("Rate of exposure to intervention") + 
  theme_minimal() + 
  theme(panel.border = element_rect("black", fill = NA, size = .5))

ggsave(fig_5, file = "results/Figure_5.pdf", width = 8, height = 4)

# Table A2 ----------------------------------------------------------------

a1 <- lm_robust(I(1-in_endline) ~ Z_common + Z_alt, data = survey, clusters = nro_cuadra)
a2 <- lm_robust(I(1-in_endline) ~ Z_common + Z_alt + as.factor(block_ID), data = survey, clusters = nro_cuadra)
a3 <- lm_robust(I(1-in_endline) ~ Z_common * Z_alt, data = survey, clusters = nro_cuadra)
a4 <- lm_robust(I(1-in_endline) ~ Z_common * Z_alt + as.factor(block_ID), data = survey, clusters = nro_cuadra)

write(texreg(list(a1, a2, a3, a4),
             include.ci = F, 
             omit.coef = c('block|(Intercept)'), 
             stars = c(.01, .05, .1), digits = 3), file = "results/Table_A2.tex")

# Table A9 ----------------------------------------------------------------

bl_trust <- group_by(survey, CUADRANTE) %>%
  summarize(mean_trust = mean(bl_trust, na.rm = T)) %>%
  mutate(mean_trustZ =(mean_trust - mean(mean_trust))/sd(mean_trust)) 
outcomes <- left_join(outcomes, bl_trust)


m1 <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ, data = outcomes, clusters = CUADRANTE)
m1a <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ + as.factor(block_ID), data = outcomes, clusters = CUADRANTE)
m1b <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ+ as.factor(block_ID) + trust_po_bl, data = outcomes, clusters = CUADRANTE)

m2 <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ, data = outcomes, clusters = CUADRANTE)
m2a <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ + as.factor(block_ID), data = outcomes, clusters = CUADRANTE)
m2b <- lm_robust(trust_beat_po ~ (Z_common + Z_alt) * mean_trustZ+ as.factor(block_ID) + trust_po_bl, data = outcomes, clusters = CUADRANTE)

m3 <- lm_robust(I(trust_beat_po - trust_police) ~ (Z_common + Z_alt) * mean_trustZ, data = outcomes, clusters = CUADRANTE)
m3a <- lm_robust(I(trust_beat_po - trust_police) ~ (Z_common + Z_alt) * mean_trustZ + as.factor(block_ID), data = outcomes, clusters = CUADRANTE)
m3b <- lm_robust(I(trust_beat_po - trust_police) ~ (Z_common + Z_alt) * mean_trustZ+ as.factor(block_ID) + trust_po_bl, data = outcomes, clusters = CUADRANTE)

m4 <- lm_robust(I(trust_police - trust_po_bl) ~ (Z_common + Z_alt) * mean_trustZ, data = outcomes, clusters = CUADRANTE)
m4a <- lm_robust(I(trust_police - trust_po_bl) ~ (Z_common + Z_alt) * mean_trustZ + as.factor(block_ID), data = outcomes, clusters = CUADRANTE)
m4b <- lm_robust(I(trust_police - trust_po_bl) ~ (Z_common + Z_alt) * mean_trustZ+ as.factor(block_ID) + trust_po_bl, data = outcomes, clusters = CUADRANTE)
library(texreg)

texreg(list(m1, m1a, m1b, m2, m2a, m2b, m3, m3a, m3b, m4, m4a, m4b), omit.coef = c("block|_bl|Intercept"), include.ci = F, digits = 3,
       custom.gof.rows = list("Block FE" = c(rep(c("", "$\\checkmark$", "$\\checkmark$"), 12/3)),
                              "Baseline outcome" = c(rep(c("", "", "$\\checkmark$"), 12/3))),
       custom.coef.names = c("Meetings", "Flyers", "Baseline beat trust (std.)", "Meetings * Beat trust (std.)", "Flyers * Beat trust (std.)"),
       custom.model.names = paste("(",as.character(c(1:12)), ")", sep = ""),
       file = "results/Table_A9.tex")


# Figure A5 ---------------------------------------------------------------

fig_a5a <- group_by(survey, nro_cuadra, Z_common, Z_alt) %>%
  summarize(attrit_rate = mean(1-in_endline)) %>%
  ggplot(aes(x = attrit_rate, col = as.factor(Z_common))) +
  stat_ecdf(lwd = .9) + 
  scale_color_manual("Meetings Arm", 
                     values = c("#a6bddb", "#1c9099"),
                     labels = c("No meetings", "Meetings")) +
  scale_x_continuous("Beat-level attrition rate") +
  scale_y_continuous("ECDF") +
  theme_minimal() +
  theme(legend.position = "bottom") +
  ggtitle("Community Policing Meeting Arm")


fig_a5b <- group_by(survey, nro_cuadra, Z_common, Z_alt) %>%
  summarize(attrit_rate = mean(1-in_endline)) %>%
  ggplot(aes(x = attrit_rate, col = as.factor(Z_alt))) +
  stat_ecdf(lwd = .9) + 
  scale_color_manual("Flyers Arm", 
                     values = c("#a6bddb", "#1c9099"),
                     labels = c("No flyers", "Flyers")) +
  scale_x_continuous("Beat-level attrition rate") +
  scale_y_continuous("ECDF") +
  theme_minimal() +
  theme(legend.position = "bottom") +
  ggtitle("Informational Flyer Arm")

ggsave(fig_a5a + fig_a5b, filename = "results/Figure_A5.pdf",
       width = 8, 
       height = 4)

# Figure A6 ---------------------------------------------------------------

key <- inner_join(select(survey, REGISTRO_BD, nro_cuadra),
                  select(outcomes, REGISTRO_BD, CUADRANTE)) %>%
  group_by(nro_cuadra, CUADRANTE) %>% summarize()


bl <- group_by(survey, nro_cuadra, Z_common) %>%
  summarize(exposure = mean(past_meeting),
            attended = mean(attended_past)) %>%
  mutate(round = "Baseline")

el <- group_by(outcomes, CUADRANTE, Z_common) %>%
  summarize(exposure = mean(any_exposure),
            attended = mean(attended)) %>%
  mutate(round = "Endline")

bl <- left_join(bl, key) %>% select(nro_cuadra, CUADRANTE, Z_common, exposure, attended, round)
el <- left_join(el, key) %>% select(nro_cuadra, CUADRANTE, Z_common, exposure, attended, round)

cuad_lev <- rbind(bl, el)

fig_a6a <- cuad_lev %>%
  group_by(Z_common, round) %>%
  mutate(mean = mean(exposure, na.rm = T)) %>%
  ggplot(aes(x = exposure, fill = as.factor(Z_common))) +
  geom_histogram(position = position_dodge(width = .02), binwidth = .05, alpha = .5) +
  facet_grid(round~.) +
  geom_vline(aes(xintercept = mean, lty = as.factor(Z_common)), lwd = 1) +
  theme_minimal() +
  scale_y_continuous("Number of beats") + 
  scale_fill_manual("Assignment", values = c("dark blue", "dark red"), labels = c("No meetings", "Meetings")) +
  scale_linetype_manual("Assignment", values = c(3,1), labels = c("No meetings", "Meetings")) +
  scale_x_continuous("Proportion of survey respondents per beat\n aware of police-community meetings") +
  theme(legend.position = "bottom")


fig_a6b <- cuad_lev %>%
  group_by(Z_common, round) %>%
  mutate(mean = mean(attended, na.rm = T)) %>%
  ggplot(aes(x = attended, fill = as.factor(Z_common))) +
  geom_histogram(position = position_dodge(width = .02), binwidth = .05, alpha = .5) +
  facet_grid(round~.) +
  geom_vline(aes(xintercept = mean, lty = as.factor(Z_common)), lwd = 1) +
  theme_minimal() +
  scale_y_continuous("Number of beats") + 
  scale_fill_manual("Assignment", values = c("dark blue", "dark red"), labels = c("No meetings", "Meetings")) +
  scale_linetype_manual("Assignment", values = c(3,1), labels = c("No meetings", "Meetings")) +
  scale_x_continuous("Proportion of survey respondents per beat\n that report attending police-community meetings") +
  theme(legend.position = "bottom")


ggsave(fig_a6a + fig_a6b, 
       file = "results/Figure_A6.pdf", 
       width = 10, height = 5)

# Figure A7 ---------------------------------------------------------------
## function generates bootstrap samples (resampling within cuadrante) and 
## estimates uptake as a function of 
btstrp_panel <- function(){
  dat <- outcomes[!is.na(outcomes$trust_po_bl),]
  dat$id <- 1:nrow(dat)
  bs <- dat[unlist(tapply(dat$id, 
                          INDEX = dat$CUADRANTE, 
                          FUN = function(x){sample(x, size = length(x), replace = T)})),]
  p1 <- as.vector(prop.table(table(bs$any_exposure, bs$trust_po_bl), margin = 1))
  p2 <- as.vector(prop.table(table(bs$attended, bs$trust_po_bl), margin = 1))
  return(c(p1, p2))
}

bs_reps <- replicate(btstrp_panel(), n = 1e3)
fig_a7 <- t(apply(bs_reps, 1, quantile, probs = c(.025, .975))) %>%
  as.data.frame() %>% 
  mutate(obs = c(as.vector(prop.table(table(outcomes$any_exposure, outcomes$trust_po_bl), margin = 1)),
                 as.vector(prop.table(table(outcomes$attended, outcomes$trust_po_bl), margin = 1))),
         outcome = rep(c("Any Exposure to Intervention", "Attended Meeting"), each = 8),
         Group = rep(c("No","Yes"), 8),
         bl_trust = rep(rep(c(1:4), each = 2),2)) %>%
  ggplot(aes(x = bl_trust, y = obs, fill = Group)) + 
  scale_fill_manual(values = c("#41b6c4", "#225ea8")) + 
  facet_grid(~outcome) + 
  geom_bar(alpha = .8,
           stat = "identity", position = position_dodge(width = .75)) + 
  ggtitle("Baseline Distribution of Trust") + 
  geom_errorbar(
    aes(ymin = `2.5%`, ymax = `97.5%`), position = position_dodge(width = .75), width = 0.05) +
  theme_minimal() + 
  scale_x_continuous("Baseline trust", breaks = 1:4) + 
  scale_y_continuous("Proportion of panel respondents") + 
  theme(legend.position = "bottom",
        panel.border = element_rect("black", fill = NA, size = .5)) 

ggsave(fig_a7, filename = "results/Figure_A7.pdf", width = 7, height = 4)


# Figure A8 ---------------------------------------------------------------

get_itts <- function(dv){
  form1 <- paste0(dv, " ~ Z_common + Z_alt")
  form2 <- paste0(form1, "+ as.factor(block_ID)")
  m1 <- lm_robust(as.formula(form1), 
                  data = survey,
                  clusters = nro_cuadra,
                  subset = !is.na(REGISTRO_BD))
  m2 <- lm_robust(as.formula(form2), 
                  data = survey,
                  clusters = nro_cuadra,
                  subset = !is.na(REGISTRO_BD))
  
  
  return(
    data.frame(
      rbind(summary(m1)$coef[2:3, 1:2],
            summary(m2)$coef[2:3, 1:2])) %>%
      mutate(estimator = rep(c("No FE", "Block FE"), each = 2),
             coef = rep(c("Meetings", "Flyers"), 2), 
             dv = dv))
}

get_cond_means <- function(dv){
  form1 <- paste0(dv, " ~ Z_common + Z_alt")
  form2 <- paste0(form1, "+ as.factor(block_ID)")
  m1 <- lm_robust(as.formula(form1), 
                  data = survey,
                  clusters = nro_cuadra,
                  subset = !is.na(REGISTRO_BD))
  
  means = c(coef(m1)[1],sum(coef(m1)[1:2]))
  ses = c(summary(m1)$coef[1,2],
          sqrt(vcov(m1)[1,1] + vcov(m1)[2,2] + 2 * vcov(m1)[1,2]))
  return(
    data.frame(means,
               ses,
               cond = c("No Meetings", "Meetings"),
               dv))
}

cms <- do.call("rbind",
               lapply(list("received_invite",
                           "heard_of",
                           "attended",
                           "any_exposure"),
                      FUN = get_cond_means))

ggsave(
  cms %>%
    mutate(outcome = plyr::mapvalues(dv,
                                     from = c("received_invite",
                                              "heard_of",
                                              "attended",
                                              "any_exposure"),
                                     to = c("Received invitation",
                                            "Heard of meetings,\n(not from invite)",
                                            "Attended meeting",
                                            "Any Exposure"))) %>%
    ggplot(aes(y = means, x = outcome, col = cond, shape = cond)) +
    geom_point(position = position_dodge(width = .5), size = 3) +
    geom_errorbar(aes(ymin = means - 1.96 * ses, ymax = means + 1.96 * ses), 
                  width = 0,
                  position = position_dodge(width = .5)) +
    geom_errorbar(aes(ymin = means - 1.64 * ses, ymax = means + 1.64 * ses), 
                  width = 0,
                  lwd = 1, 
                  position = position_dodge(width = .5)) +
    theme_minimal() +
    scale_x_discrete("Outcome Measure") +
    scale_y_continuous("Condition mean", limits = c(0, .45)) +
    scale_color_manual("Treatment Condition", values = c("#a6bddb", "#1c9099")) +
    scale_shape_manual("Treatment Condition", values = c(15, 16)) +
    theme(legend.position = "bottom") +
    ggtitle("Exposure to Treatment, Panel Sample") + 
    coord_flip(),
  file = "results/Figure_A8.pdf", 
  width = 5, height = 5
)


# Figure A9 ---------------------------------------------------------------
## Per the glmnet documentation, the results of lasso are random due to the use of cross validation to select lambda. These results change slightly even with the seed on line 14. 
## Our principal finding: the ranking (#2) and approximate magnitude of "highest trust in police" as a predictor of either measure of participation does not change substantially.

bl_vars <- select(survey, REGISTRO_BD, male:bl_trust,race, COMUNA_LOCALIDAD, Z_common, Z_alt, nro_cuadra)
outcome_vars <- select(outcomes, REGISTRO_BD, any_exposure, attended)
lasso_vars <- left_join(outcome_vars, bl_vars) %>%
  filter(!is.na(REGISTRO_BD) & !is.na(bl_trust))

mod.mat <- model.matrix(~-1+ male + as.factor(estrato) + as.factor(estrato_servicios) + 
                          as.factor(education) + migrant + as.factor(decade) + race + COMUNA_LOCALIDAD + 
                          pleb + not_partisan + centro_democratico + farc + would_vote + pres_vote_intention + as.factor(bl_trust) + 
                          Z_common + Z_alt, data = lasso_vars)
X <- scale(mod.mat)
X <- makeX(as.data.frame(X), na.impute = T)
any.exposure <- lasso_vars$any_exposure
attended <- lasso_vars$attended


lasso1a <- glmnet(X, any.exposure, alpha = 1, intercept = F)
lasso_cv1a <- cv.glmnet(X, any.exposure, alpha = 1, intercept = F, k = 10)
coefs_ae <- coef(lasso1a, s = lasso_cv1a$lambda.min)

lasso2a <- glmnet(X, attended, alpha = 1, intercept = F)
lasso_cv2a <- cv.glmnet(X, attended, alpha = 1, intercept = F, k = 10)
coefs_att <- coef(lasso2a, s = lasso_cv2a$lambda.min)



labs <- read.csv("data/lasso_var_key.csv")
fig_a9a <- left_join(data.frame(factor = coefs_ae@Dimnames[[1]][coefs_ae@i+1], 
                      est = coefs_ae@x), labs) %>%
  mutate(col = ifelse(label %in% c("Meetings treatment", "Flyers treatment"), "Treatment", ifelse(grepl(label, pattern = "trust"), "Baseline trust", "Other")),
         col = factor(col, levels = c("Baseline trust", "Treatment", "Other"))) %>%
  ggplot(aes(x = est, y = reorder(label, est), col = col, shape = col)) + 
  geom_point(size = 2) + 
  scale_x_continuous("LASSO coefficient estimate on 'Any Exposure'") + 
  scale_y_discrete("") + 
  theme_minimal() + 
  scale_color_manual("Variable category", values = c("dark red", "dark blue", "gray")) +
  scale_shape_manual("Variable category", values =c(17, 16, 15)) + 
  geom_vline(xintercept = 0, lty = 3, col = "red") + 
  theme(legend.position = "bottom")



fig_a9b <- left_join(data.frame(factor = coefs_att@Dimnames[[1]][coefs_att@i+1], 
                      est = coefs_att@x), labs) %>%
  mutate(col = ifelse(label %in% c("Meetings treatment", "Flyers treatment"), "Treatment", ifelse(grepl(label, pattern = "trust"), "Baseline trust", "Other")),
         col = factor(col, levels = c("Baseline trust", "Treatment", "Other"))) %>%
  ggplot(aes(x = est, y = reorder(label, est), col = col, shape = col)) + 
  geom_point(size = 2) + 
  scale_x_continuous("LASSO coefficient estimate on 'Attended Meeting'") + 
  scale_y_discrete("") + 
  theme_minimal() + 
  scale_color_manual("Variable category", values = c("dark red", "dark blue", "gray")) +
  scale_shape_manual("Variable category", values =c(17, 16, 15)) + 
  geom_vline(xintercept = 0, lty = 3, col = "red") + 
  theme(legend.position = "bottom")




ggsave(fig_a9a + fig_a9b + plot_layout(nrow = 2), file = "results/Figure_A9.pdf",
width = 9, height = 7)




# Figure A10 --------------------------------------------------------------
lasso_vars2 <- lasso_vars %>% mutate(trust_rescale = plyr::mapvalues(bl_trust, 1:4, 0:3/3),
                                     Z = 1.5 * Z_common + .75 * Z_alt) # creates four arbitrary treatment indicators

get_trust_coefs <- function(mod){
  return(summary(mod)$coef[grep(names(coef(mod)), pattern = "trust"), 1:2])
}

m1 <- lm_robust(any_exposure ~ trust_rescale, data = lasso_vars2, cluster = nro_cuadra)
m2 <- lm_robust(any_exposure ~ as.factor(bl_trust), data = lasso_vars2, cluster = nro_cuadra)
m3 <- lm_robust(attended ~ trust_rescale, data = lasso_vars2, cluster = nro_cuadra)
m4 <- lm_robust(attended ~ as.factor(bl_trust), data = lasso_vars2, cluster = nro_cuadra)

univariate <- do.call("rbind", lapply(list(m1, m2, m3, m4), get_trust_coefs)) %>%
  as.data.frame() %>%
  mutate(Trust = rep(c("0-1 scale", 
                       "Trust a little", "Somewhat trust", "Trust a lot"), 2),
         Panel = rep(c("Linearized", rep("Dummies", 3)),2),
         Outcome = rep(c("Any exposure", "Attended"), each = 4),
         var = "Univariate")


get_ests <- function(var, var_name){
  form1 = paste0("any_exposure ~ trust_rescale + as.factor(", var, ")")
  form2 = paste0("any_exposure ~ as.factor(bl_trust) + as.factor(", var, ")")
  form3 = paste0("attended ~ trust_rescale + as.factor(", var, ")")
  form4 = paste0("attended ~ as.factor(bl_trust) + as.factor(", var, ")")
  
  m1 <- lm_robust(as.formula(form1), data = lasso_vars2, cluster = nro_cuadra)
  m2 <- lm_robust(as.formula(form2), data = lasso_vars2, cluster = nro_cuadra)
  m3 <- lm_robust(as.formula(form3), data = lasso_vars2, cluster = nro_cuadra)
  m4 <- lm_robust(as.formula(form4), data = lasso_vars2, cluster = nro_cuadra)
  
  return(do.call("rbind", lapply(list(m1, m2, m3, m4), get_trust_coefs)) %>%
           as.data.frame() %>%
           mutate(Trust = rep(c("0-1 scale", 
                                "Trust a little", "Somewhat trust", "Trust a lot"), 2),
                  Panel = rep(c("Linearized", rep("Dummies", 3)),2),
                  Outcome = rep(c("Any exposure", "Attended"), each = 4),
                  var = var_name))
}

cov.adj <- do.call("rbind", list(
  get_ests("male", "Gender"),
  get_ests("estrato", "Stratum"),
  get_ests("education", "Education"),
  get_ests("migrant", "Migrant"),
  get_ests("decade", "Age (decade bins)"),
  get_ests("race", "Race"),
  get_ests("not_partisan", "Non-partisan indicator"),
  get_ests("centro_democratico", "Centro Democratico partisan"),
  get_ests("farc", "FARC partisan"),
  get_ests("would_vote", "Would vote in upcoming\npresidential election"),
  get_ests("pres_vote_intention", "Presidential Vote Intention"),
  get_ests("Z", "Treatments")))


levs <- c(sort(unique(cov.adj$var)), "Univariate")
cols <- c(brewer.pal(12, "Set3"), "black")
fig_a10 <- rbind(univariate, cov.adj) %>%
  mutate(Trust = factor(Trust, levels = unique(Trust)),
         `Model (covariate)` = factor(var, levels = c(levs))) %>%
  ggplot(aes(y = Trust, x = Estimate, col = `Model (covariate)`)) + 
  geom_point(position = position_dodge(width =.75)) + 
  facet_grid(Panel ~ Outcome, space = "free", scales = "free") +
  geom_errorbarh(aes(xmin = Estimate - 1.96 * `Std. Error`,
                     xmax = Estimate + 1.96 * `Std. Error`), height = 0,
                 position = position_dodge(width = .75)) + 
  scale_color_manual(values = cols) + 
  geom_vline(xintercept = 0, lty = 3) + 
  theme_minimal() + 
  theme(  panel.border = element_rect(colour = "black", fill=NA, size=1))

ggsave(fig_a10, file = "results/Figure_A10.pdf", width = 9, height = 5)

# Figure A11 --------------------------------------------------------------

fig_a11a <- ggplot(outcomes, aes(x = age, y = trust_po_bl)) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F) + 
  scale_x_continuous("Age at baseline") + 
  scale_y_continuous("Prior: Baseline trust in police") + 
  theme_minimal()

fig_a11b <- filter(outcomes, Z_common == 1) %>% 
  ggplot(aes(x = age, y = any_exposure)) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F) + 
  scale_x_continuous("Age at baseline") + 
  scale_y_continuous("Any exposure to community policing\n(treatment group)") + 
  theme_minimal()


fig_a11c <- ggplot(outcomes, aes(x = age, y = trust_police-trust_po_bl, col = as.factor(Z_common))) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F)+ 
  scale_x_continuous("Age at baseline") + 
  scale_y_continuous("Change in trust in police (endline-baseline)") + 
  scale_color_manual("Treatment group", 
                     values = c("#a6bddb", "#2b8cbe"), 
                     labels = c("No meetings", "Meetings")) +
  theme_minimal() + 
  theme(legend.position = "bottom")

fig_a11d <- ggplot(outcomes, aes(x = migrant, y = trust_po_bl)) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F, method = "lm") + 
  scale_x_continuous("Internal migrant") + 
  scale_y_continuous("Prior: Baseline trust in police") + 
  theme_minimal()

fig_a11e <- filter(outcomes, Z_common == 1) %>% 
  ggplot(aes(x = migrant, y = any_exposure)) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F, method = "lm") + 
  scale_x_continuous("Internal migrant") + 
  scale_y_continuous("Any exposure to community policing\n(treatment group)") + 
  theme_minimal()

fig_a11f <- ggplot(outcomes, aes(x = migrant, y = trust_police-trust_po_bl, col = as.factor(Z_common))) + 
  geom_jitter(size = .2, height = .2) + 
  geom_smooth(se = F, method = "lm")+ 
  scale_x_continuous("Internal migrant") + 
  scale_y_continuous("Change in trust in police (endline-baseline)") + 
  scale_color_manual("Treatment group", 
                     values = c("#a6bddb", "#2b8cbe"), 
                     labels = c("No meetings", "Meetings")) +
  theme_minimal() +
  theme(legend.position = "bottom")

library(patchwork)
ggsave(fig_a11a + fig_a11b + fig_a11c + fig_a11d + fig_a11e + fig_a11f + plot_layout(ncol = 3),
       file = "results/Figure_A11.pdf", width = 11, height = 8)

# Figure A12 --------------------------------------------------------------


fig_a12 <- data.frame(po_index_bl = rep(outcomes$po_index_bl, 2),
                 d = c(outcomes$any_exposure, outcomes$attended),
                 outcome = rep(c("Any exposure", "Attended Meeting"), each = nrow(outcomes))) %>%
  ggplot(aes(x = po_index_bl, y = d)) +
  geom_jitter(height = .02) +
  geom_smooth() + 
  facet_grid(~outcome) + 
  xlab("Baseline police quality index") + 
  ylab("Proportion of respondents in panel") + 
  theme_minimal()

ggsave(fig_a12, file = "results/Figure_A12.pdf", width = 8, height = 4)

# Figure A15 --------------------------------------------------------------
# subset to respondents in panel
panel <- filter(outcomes, !is.na(trust_po_bl))

# estimates extreme value bounds on itt
evb <- c(4-mean(panel$trust_po_bl[panel$Z_common ==0]),
         1-mean(panel$trust_po_bl[panel$Z_common ==0]))

names(evb) <- c("ub", "lb")

# Rate of compliance
comp.ae <- lm_robust(any_exposure ~ -1 + as.factor(Z_common), data= panel)$coef
comp.am <- lm_robust(attended ~ -1 + as.factor(Z_common), data= panel)$coef
bl.average <- tapply(panel$trust_po_bl, panel$Z_common, FUN = mean)

evb.comp1 <- c(comp.ae[2] * 4 + (1-comp.ae[2]) * bl.average[2] -
                 (comp.ae[1] * 4 + (1-comp.ae[1]) * bl.average[1]),
               comp.ae[2] * 1 + (1-comp.ae[2]) * bl.average[2] -
                 (comp.ae[1] * 1 + (1-comp.ae[1]) * bl.average[1]))


evb.comp2 <- c(comp.am[2] * 4 + (1-comp.am[2]) * bl.average[2] -
                 (comp.am[1] * 4 + (1-comp.am[1]) * bl.average[1]),
               comp.am[2] * 1 + (1-comp.am[2]) * bl.average[2] -
                 (comp.am[1] * 1 + (1-comp.am[1]) * bl.average[1]))

panel %<>% mutate(ub.ae = ifelse(any_exposure == 1, 4, trust_po_bl),
                  lb.ae = ifelse(any_exposure == 1, 1, trust_po_bl),
                  ub.at = ifelse(attended == 1, 4, trust_po_bl),
                  lb.at = ifelse(attended == 1, 1, trust_po_bl))

evb.comp1.s <- c(lm_robust(ub.ae ~ Z_common, data = panel)$coef[2],
                 lm_robust(lb.ae ~ Z_common, data = panel)$coef[2])

evb.comp2.s <- c(lm_robust(ub.at ~ Z_common, data = panel)$coef[2],
                 lm_robust(lb.at ~ Z_common, data = panel)$coef[2])


x <- seq(from = 0, to = .55, by = .01)
fig_a15 <- data.frame(ub = c(x * evb.comp1[1]/diff(comp.ae),
                        x * evb.comp1.s[1]/diff(comp.ae),
                        x * evb.comp2[1]/diff(comp.am),
                        x * evb.comp2.s[1]/diff(comp.am)),
                 lb = c(x * evb.comp1[2]/diff(comp.ae),
                        x * evb.comp1.s[2]/diff(comp.ae),
                        x * evb.comp2[2]/diff(comp.am),
                        x * evb.comp2.s[2]/diff(comp.am)),
                 compliance = rep(x, 4), 
                 plot = rep(c("Any Exposure", "Attended Meeting"), each = 112),
                 selection = rep(c("Random", "Observed", "Random", "Observed"), each = 56)) %>%
  mutate(observed = ifelse(plot == "Any Exposure", diff(comp.ae), diff(comp.am))) %>%
  ggplot(aes(x = compliance, ymin = lb, ymax = ub, group = selection, fill = selection)) + 
  geom_ribbon(alpha = 0.4) + 
  facet_grid(~plot) +
  geom_vline(aes(xintercept = observed), lty = 2) + 
  theme_minimal() + 
  scale_fill_manual("Selection", values = c("#fdb863", "#b2abd2")) + 
  xlab(expression(paste(ITT[D], ": Effect on compliance"))) + 
  ylab("Bounds on the ITT") + 
  theme(legend.position = "bottom")

ggsave(fig_a15, file = "results/Figure_A15.pdf", width = 7, height = 4)

# Figure A16 --------------------------------------------------------------

# function used to construct index

Z_score <- function(var, Z) {(var - mean(var[Z==0], na.rm = T))/sd(var[Z==0], na.rm = T)}

greedy_indexing <- function(vars, Z){
  non_missing = rowSums(!is.na(vars))
  Z_scores = rowSums(vars, na.rm = T)
  index_Z = Z_score(Z_scores/non_missing, Z)
  return(index_Z)
}

# impute higher outcomes for relevant individuals under each scenario
outcomes %<>% 
  mutate(tc1 = ifelse(trust_po_bl  == 4 & trust_police == 4 & Z_common == 1, 5, trust_police), 
         tc2 = ifelse(trust_po_bl  == 4 & trust_police == 4 & any_exposure == 1, 5, trust_police),
         tc3 = ifelse(trust_po_bl  == 4 & trust_police == 4 & attended == 1, 5, trust_police))

uncensor_index <- function(vector){
  aa <- ifelse(vector == 1 & outcomes$act_appropriate_po == 5 & outcomes$act_appropriate_po_bl == 5,
               6, outcomes$act_appropriate_po)
  as <- ifelse(vector == 1 & outcomes$act_seriously_po == 5 & outcomes$act_seriously_po_bl == 5,
               6, outcomes$act_seriously_po)
  es <- ifelse(vector == 1 & outcomes$equal_service_po == 5 & outcomes$equal_service_po_bl == 5,
               6, outcomes$equal_service_po)
  co <- ifelse(vector == 1 & outcomes$corrupt_po == 5 & outcomes$corrupt_po_bl == 5,
               6, outcomes$corrupt_po)
  ci <- ifelse(vector == 1 & outcomes$capacity_investigate_po_bl == 5 & outcomes$capacity_investigate_po == 5,
               6, outcomes$capacity_investigate_po)
  cr <- ifelse(vector == 1 & outcomes$capacity_respond_po_bl == 5 & outcomes$capacity_respond_po == 5,
               6, outcomes$capacity_respond_po)
  cens_po_index <- greedy_indexing(cbind(aa, as, es, co, ci, cr), 1 *(outcomes$Z_common + outcomes$Z_alt > 0))
  return(cens_po_index)
}

outcomes %<>% mutate(
  pc1 = uncensor_index(Z_common),
  pc2 = uncensor_index(any_exposure),
  pc3 = uncensor_index(attended)
)

# Estimate ITTs with each imputed outcome for those in panel sample
baseline <- lm_robust(trust_police ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                      data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound1 <- lm_robust(tc1 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound2 <- lm_robust(tc2 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound3 <- lm_robust(tc3 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))

baseline2 <- lm_robust(po_index ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                       data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound4 <- lm_robust(pc1 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound5 <- lm_robust(pc2 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))
bound6 <- lm_robust(pc3 ~ Z_common + Z_alt + as.factor(block_ID) + trust_po_bl, 
                    data = outcomes, clusters = CUADRANTE, subset = !is.na(REGISTRO_BD))

ggsave(do.call("rbind",
               lapply(list(baseline, bound1, bound2, bound3,
                           baseline2, bound4, bound5, bound6), FUN = function(x){summary(x)$coef[2,1:2]})) %>%
         as.data.frame() %>%
         mutate(estimate = rep(c("Observed", "Censoring for all of\ntreatment group", "Censoring for anyone\nexposed to meetings",
                                 "Censoring for all\nmeeting attendees"),2),
                estimate = factor(estimate, rev(unique(estimate))),
                outcome = rep(c("Trust in police", "Police quality index"), each = 4)) %>%
         ggplot(aes(x = Estimate, y= estimate)) + 
         xlab("ITT estimate") + 
         ylab("Scenario") + 
         geom_point(size = 2) + 
         facet_grid(~outcome) + 
         geom_errorbar(aes(xmin = Estimate - 1.96 * `Std. Error`, xmax = Estimate + 1.96 * `Std. Error`), width = 0.1) + 
         geom_vline(xintercept = 0, lty = 3) + 
         theme_minimal(), width = 9, height = 3, file = "results/Figure_A16.pdf")



# Figure A17 --------------------------------------------------------------

# Examine individuals in treatment cuadrantes
sset <- filter(panel, Z_common == 1) %>%
  mutate(id = 1:n())

sset %<>% mutate(id = 1:n())


options(na.action = "na.omit")

# predict outcomes for individuals 
yhat <- lm(trust_police ~ -1 + as.factor(trust_po_bl), data =sset, subset = any_exposure == 0)
df_at <- model.matrix(~-1 + as.factor(trust_po_bl), data = sset[sset$attended == 1,]) 
df_ae <- model.matrix(~-1 + as.factor(trust_po_bl), data = sset[sset$attended == 0 & sset$any_exposure == 1,]) 

# Neither
mean(predict(yhat))
#Exposure
mean(df_ae %*% coef(yhat))
# Attendees
mean(df_at %*% coef(yhat))

bootstrap <- function(){
  bs.dat <- sset[unlist(tapply(sset$id, INDEX = sset$CUADRANTE, FUN = function(x){sample(x, size = length(x), replace = T)})),]
  
  yhat <- lm(trust_police ~ -1 + as.factor(trust_po_bl), data =bs.dat, subset = any_exposure == 0)
  df_at <- model.matrix(~-1 + as.factor(trust_po_bl), data = bs.dat[bs.dat$attended == 1,]) 
  df_ae <- model.matrix(~-1 + as.factor(trust_po_bl), data = bs.dat[bs.dat$attended == 0 & bs.dat$any_exposure == 1,]) 
  
  # Neither
  out <- c(mean(predict(yhat)), mean(df_ae %*% coef(yhat)), mean(df_at %*% coef(yhat)))
  names(out) <- c("neither", "heard of", "attended")
  return(out)
}

sims <- replicate(n = 1e3, expr = bootstrap())

ggsave(data.frame(sims = as.vector(sims), 
                  group = rep(c("Neither", "Heard of Meetings", "Attended Meeting"), 1e3)) %>%
         mutate(obs = ifelse(group == "Neither", mean(sset$trust_police[sset$any_exposure == 0], na.rm = T),
                             ifelse(group == "Heard of Meetings", mean(sset$trust_police[sset$attended == 0 & sset$any_exposure == 1], na.rm = T),
                                    mean(sset$trust_police[sset$attended == 1], na.rm = T)))) %>%
         mutate(group = factor(group, levels = unique(group))) %>%
         filter(group != "Neither") %>%
         group_by(group) %>%
         mutate(p = paste0("p = ", round(2 * mean(obs < sims), 3))) %>%
         ggplot(aes(x = sims, label = p)) + 
         facet_grid(~group) + 
         geom_histogram() + 
         geom_text(aes(x= obs + .08), y = 140, col = "red") + 
         geom_vline(aes(xintercept = obs), col = "red", lwd = .8) + 
         ggtitle("Bootstrapped predictive distribution of endline trust in police assuming\nmean regression at rate of mean regression among unexposed citizens") + 
         scale_x_continuous("Average trust in police") + 
         theme_minimal(), file = "results/Figure_A17.pdf", width = 7, height = 4.25)




